Report Gruppe E

Introduction and data

Project Motivation

1999 wurde in der Forschung (Collins et al.) ermittelt, dass die ASA-Klassifizierung des körperlichen Zustands, das Alter, Bluttransfusionen während der Operation und die Dauer von Operationen einen Aufschluss darüber geben, ob der Krankenhausaufenthalt des Patienten länger als gewöhnlich ausfallen wird. Des Weiteren gibt es zahlreiche Studien, bspw. die Forschung von Biber et al. (2012), Singler et al. (2013) oder Ogliari et al. (2022) die festgestellt haben, dass es eine Korrelation zwischen dem Alter von Patienten und der Verweildauer in Notaufnahmen gab.

Dieser Frage möchten wir gerne nachgehen und versuchen, die Aufenthaltsdauer im Krankenhaus vorherzusagen.

Data

Uns liegen Daten von nicht-kardiochirurgischen Patienten vor, welche sich von August 2016 bis Juni 2017 einer Routine- oder Notoperation am Seoul National University Hospital, Seoul, Korea, unterzogen. Von den 7.051 in Frage kommenden Fällen wurden 663 qualitativ unzureichende Fälle ausgeschlossen. Schließlich wurden 6.388 Fälle (91 %) in den Datensatz aufgenommen. Die Daten wurden nicht vorverarbeitet, da das reale Rauschen in den Vitaldaten für die Entwicklung praktischer Überwachungsalgorithmen sehr wichtig ist.

Research Question

Wir stellen uns die Frage, welche dem Krankenhaus vorliegenden Daten die Aufenthaltsdauer von Patienten beeinflusst. Da der Datensatz des Seoul National University Hospital (Seoul, Korea) eine Vielzahl von möglichen Variablen enthält werden wir uns vorerst mit den aus der Literatur erwähnten Zusammenhängen beschäftigen.

Code

import pandas as pd
import altair as alt
import numpy as np
import joblib
import seaborn as sns
import matplotlib.pyplot as plt

from sklearn.feature_selection import SelectFromModel

from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.inspection import permutation_importance

from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV
from sklearn.linear_model import Lasso

from statsmodels.formula.api import ols

from sklearn import datasets, ensemble
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
alt.data_transformers.disable_max_rows()

DataTransformerRegistry.enable('default')

Data preparation

Daten importieren

Code

df = pd.read_csv("https://raw.githubusercontent.com/Maarv99/Project-applied-statistics/main/data/raw/clinical_data.csv")
df_orig = df.copy()

Wir berechnen eine weitere Variable “op_duartion”. Sie zeigt die Dauer einer Operation an. Dafür subtrahieren wir von der Variablen “opend” die Variable “opstart”.

Code

df['op_duration']=df['opend']-df['opstart']

Wir benennen einige Variablen um:

• “caseend” bezeichnet das Ende eines Falls, vom Startpunkt 0, daher bennen wir diese Variable “record_duration”
• Die Variable “icu_days” bezeichnet die Aufenthaltsdauer im Krankenhaus (icu), daher bennen wir diese Variable “length_of_stay”
• Da wir uns unter der Variable “intraop_ebl” wenig vorstellen konnten, benennen wir diese in “estimated_blood_loss” um.

Code

df=df.rename(columns={'caseend':'recording_duration','icu_days':'length_of_stay','intraop_ebl':'estimated_blood_loss'})

Wir ändern nun den Datentyp der Variable “estimated_blood_loss”, sodass diese korrekt erkannt wird.

Code

df['estimated_blood_loss'] = df['estimated_blood_loss'].fillna(0).astype(int)

Die Spalte “age” enthält u.a. die Zeichen “>89” und “0.7”. Daher wird das Größer-als-Zeichen entfernt und die Werte gerundet, damit es zu int umwandelt werden kann.

Code

df['age'] = df['age'].str.replace('>89', '89').astype(float).round().astype(int)

Wir prüfen unsere Daten auf Zellen ohne Werte.

Code

print("Missing values in 'estimated_blood_loss':",df['estimated_blood_loss'].isnull().sum())
print("Missing values in 'op_duration':",df['op_duration'].isnull().sum())
print("Missing values in 'length_of_stay':",df['length_of_stay'].isnull().sum())
print("Missing values in 'bmi':",df['bmi'].isnull().sum())
print("Missing values in 'asa':",df['asa'].isnull().sum())
print("Missing values in 'age':",df['age'].isnull().sum())
print("Missing values in 'recording_duration':",df['recording_duration'].isnull().sum())
print("Missing values in 'sex':",df['sex'].isnull().sum())

Missing values in 'estimated_blood_loss': 0
Missing values in 'op_duration': 0
Missing values in 'length_of_stay': 0
Missing values in 'bmi': 0
Missing values in 'asa': 133
Missing values in 'age': 0
Missing values in 'recording_duration': 0
Missing values in 'sex': 0

Lediglich die Spalte ‘asa’ hat 133 NAs. Wir werden diese Variable mit Hilfe eines Boxplots genauer betrachten um zu entscheiden wie wir damit umgehen.

Code

# Boxplot 
box = alt.Chart(df).mark_boxplot().encode(
    x='asa',
)

box

Wir werden die NAs durch den Median ersetzen und prüfen, ob dies erfolgreich war.

Code

df['asa'] = df['asa'].fillna(df['asa'].median())
print("Missing values in 'asa':",df['asa'].isnull().sum())

Missing values in 'asa': 0

Wir ändern kategorische Werte in numerische und lassen uns die Daten einmal anzeigen:

Code

df = df[['estimated_blood_loss', 'op_duration', 'length_of_stay', 'bmi', 'asa', 'age', 'recording_duration', 'sex']]
df = pd.get_dummies(df)
df

	estimated_blood_loss	op_duration	length_of_stay	bmi	asa	age	recording_duration	sex_F	sex_M
0	0	8700	0	26.3	2.0	77	11542	0	1
1	50	12900	0	19.6	2.0	54	15741	0	1
2	0	1920	0	24.4	1.0	62	4394	0	1
3	0	15300	1	20.5	2.0	74	20990	0	1
4	2600	17700	13	20.4	3.0	66	21531	0	1
...	...	...	...	...	...	...	...	...	...
6383	100	12000	0	24.2	1.0	64	15248	0	1
6384	100	17100	0	24.6	2.0	69	20643	0	1
6385	100	14700	0	18.8	1.0	61	19451	1	0
6386	0	9300	0	22.9	1.0	24	12025	1	0
6387	0	6900	0	22.9	2.0	47	10249	1	0

6388 rows × 9 columns

Wir schauen uns noch die Verteilung der Daten an:

Code

df=df[['recording_duration','age','asa','bmi','length_of_stay','op_duration','estimated_blood_loss','sex_F','sex_M']]


y_label = 'length_of_stay'

features = ['recording_duration',
            'age',
            'asa',
            'bmi',
            'op_duration',
            'estimated_blood_loss', 
            'sex_F']

X = df[features]
y = df[y_label]

Code

alt.Chart(df).mark_circle().encode(
    alt.X(alt.repeat("column"), type='quantitative'),
    alt.Y(alt.repeat("row"), type='quantitative')
).properties(
    width=150,
    height=150
).repeat(
    row=['length_of_stay'],
    column=['age','asa','bmi','op_duration','estimated_blood_loss','sex_F','sex_M']
).interactive()

Von den 6.388 Werten lassen sich zwei Ausreißer erkennen. Im folgenden Boxplot sind diese noch einmal deutlich zu erkennen. Um unser Modell nicht von diesen zu beeinflussen, werden wir diese Ausreißer heraus nehmen. Wir nehmen Aufzeichnung mit einer “length_of_stay” größer als 50 Tagen aus den Daten heraus.

Code

# Boxplot 
box = alt.Chart(df).mark_boxplot().encode(
    x='length_of_stay',
)

box

Code

df = df[df['length_of_stay']<50]

Variablen

Aufgrund unserer Projektmotivation heraus nutzen wir die dort erwähnten Variablen. Zusätzlich nehmen wir noch die Variable “bmi” hinzu. Diese wird zwar nicht erwähnt, zeigt jedoch auf eine andere Art als die ASA-Klassifizierung (Variable “asa”) auch den Körperlichen Zustand eines Patienten. Aus der Projektmotivation ist auch zu entnehmen, dass wir die Krankenhausaufenhaltsdauer anhand dieser Prädiktoren bestimmen möchten.

Code

df=df[['recording_duration','age','asa','bmi','length_of_stay','op_duration','estimated_blood_loss','sex_F','sex_M']]


y_label = 'length_of_stay'

features = ['recording_duration',
            'age',
            'asa',
            'bmi',
            'op_duration',
            'estimated_blood_loss', 
            'sex_F']

X = df[features]
y = df[y_label]

Wir möchten uns nun die Daten einmal einzeln anschauen:

Code

alt.Chart(df).mark_bar().encode(
    alt.X(alt.repeat("column"), type="quantitative", bin=True),
    y='count()',
).properties(
    width=150,
    height=150
).repeat(
    column=['length_of_stay','age','asa','bmi','op_duration','estimated_blood_loss','sex_F','sex_M']
)

Nun Betrachten wir die Korrelation der Variablen:

Code

df[['length_of_stay','age','asa','bmi','op_duration','estimated_blood_loss','sex_F','sex_M']].corr()
corr = df.corr()
corr[y_label].sort_values(ascending=False)
corr.style.background_gradient(cmap='Blues')

	recording_duration	age	asa	bmi	length_of_stay	op_duration	estimated_blood_loss	sex_F	sex_M
recording_duration	1.000000	0.057682	0.127184	-0.088780	0.213018	0.989792	0.305531	-0.147491	0.147491
age	0.057682	1.000000	0.210617	0.067234	-0.037458	0.051759	0.019490	-0.128332	0.128332
asa	0.127184	0.210617	1.000000	-0.069810	0.296588	0.118179	0.141434	-0.152637	0.152637
bmi	-0.088780	0.067234	-0.069810	1.000000	-0.091738	-0.087440	-0.056082	0.011894	-0.011894
length_of_stay	0.213018	-0.037458	0.296588	-0.091738	1.000000	0.212291	0.270221	-0.039277	0.039277
op_duration	0.989792	0.051759	0.118179	-0.087440	0.212291	1.000000	0.314199	-0.142918	0.142918
estimated_blood_loss	0.305531	0.019490	0.141434	-0.056082	0.270221	0.314199	1.000000	-0.070212	0.070212
sex_F	-0.147491	-0.128332	-0.152637	0.011894	-0.039277	-0.142918	-0.070212	1.000000	-1.000000
sex_M	0.147491	0.128332	0.152637	-0.011894	0.039277	0.142918	0.070212	-1.000000	1.000000

Wir können eine starke Interkorrelation zwischen “op-duration” und “recording_duration” erkennen. Wir erkennen eine logisch nachvollziehbare, eindeutige, negative Übereinstimmung zwischen den Variablen “sex_F” und “sex_M”. Dass diese Variablen unser Modell nicht stören, werden wir uns für eine entscheiden und nur die Variable “sex_F” in das Modell einbeziehen und definieren daher X und Y neu. Die Korrelation der Prädiktoren mit der Variablen “length_of_stay” ist sehr gering.

Code

df=df[['recording_duration','age','asa','bmi','length_of_stay','op_duration','estimated_blood_loss','sex_F']]

Nun teilen wir den Datensatz in einen Trainings- und einen Testdatensatz auf.

Code

X_train, X_test, y_train, y_test = train_test_split(X, y, 
                                                    test_size=0.3,
                                                    random_state=42)

Methodology

Lineare Regression

Wir beginnen nun ein Regressionsmodell aufzustellen

Code

reg = LinearRegression()

Nun erstellen wir die Kreuzvalidierung der Daten mit 5 Folds

Code

scores = cross_val_score(reg, X_train, y_train, cv=5, scoring='neg_mean_squared_error') *-1
# store cross-validation scores
df_scores = pd.DataFrame({"lr": scores})

# reset index to match the number of folds
df_scores.index += 1

alt.Chart(df_scores.reset_index()).mark_line(
     point=alt.OverlayMarkDef()
).encode(
    x=alt.X("index", bin=False, title="Fold", axis=alt.Axis(tickCount=5)),
    y=alt.Y("lr", aggregate="mean", title="Mean squared error (MSE)")
)

Im nächsten Schritt passen wir das Modell an die Trainingsdaten an.

Code

reg.fit(X_train, y_train)

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Wir lassen uns nun die Regressionsgerade unsereres Modells berechnen.

Code

intercept = pd.DataFrame({
    "Name": ["Intercept"],
    "Coefficient":[reg.intercept_]}
    )

# make a slope table
slope = pd.DataFrame({
    "Name": features,
    "Coefficient": reg.coef_}
)

# combine estimates of intercept and slopes
table = pd.concat([intercept, slope], ignore_index=True, sort=False)

round(table, 3)

	Name	Coefficient
0	Intercept	-0.175
1	recording_duration	0.000
2	age	-0.016
3	asa	0.980
4	bmi	-0.035
5	op_duration	0.000
6	estimated_blood_loss	0.001
7	sex_F	0.131

Nun möchten wir sehen, wie gut unser Modell ist. Dafür prognostizieren wir Werte für die Aufenthaltsdauer (“length_of_stay”) auf Basis der Testdaten. Um das Modell bewerten zu können berechnen wir Gütemaße für die prognostizierten Werte. Diese sind R-squared (R²), Adjusted R-Squared (Adjusted R²), den Rooted-Mean-Squared-Error (RMSE), den Mean-Squared-Error (MSE) und den Mean-Absolute-Error (MAE).

Code

y_pred = reg.predict(X_test)
print("R²:",r2_score(y_test, y_pred).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test, y_pred))*((len(X_test)-1)/(len(X_test)-len(X_test.columns)-1))).round(4))
print("RMSE:", mean_squared_error(y_test, y_pred, squared=False).round(3))
print("MSE:",mean_squared_error(y_test, y_pred).round(3))
print("MAE:",mean_absolute_error(y_test, y_pred).round(3))

R²: 0.115
Adjusted R²: 0.1117
RMSE: 1.944
MSE: 3.777
MAE: 0.803

Nun möchten wir prüfen, ob unser Modell mit weniger Variablen besser wird, da beispielsweise eine starke Interkorrelation zwischen op-duration und recording_duration erkennbar ist. Wir gehen dafür entsprechend der Vorwärtseliminierung vor. Daher betrachten wir zuerst, wie wichtig welche Variable für unser Modell ist. Wir werden im Anschluss unser Modell der Reihe nach mit den wichtigsten Variablen aufbauen und die Gütemaße berechnen, bis diese nicht mehr besser werden. Hauptausschlaggebend wird für uns das adjusted R² sein.

Code

importance = np.abs(reg.coef_)

df_imp = pd.DataFrame({"coeff": importance, 
                       "name": features})
alt.Chart(df_imp).mark_bar().encode(
    x="coeff",
    y=alt.Y("name", sort='-x')
)

Code

y_label2 = 'length_of_stay'

features2 = [#'recording_duration',
            #'age',
            'asa',
            #'bmi',
            #'op_duration',
            #'estimated_blood_loss', 
            #'sex_F',
            ]

X2 = df[features2]
y2 = df[y_label2]

X_train2, X_test2, y_train2, y_test2 = train_test_split(X2, y2, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train2, y_train2)
y_pred2 = reg.predict(X_test2)

print("R²:",r2_score(y_test2, y_pred2).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test2, y_pred2))*((len(X_test2)-1)/(len(X_test2)-len(X_test2.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test2, y_pred2, squared=False).round(3))
print("MSE:",mean_squared_error(y_test2, y_pred2).round(3))
print("MAE:",mean_absolute_error(y_test2, y_pred2).round(3))

R²: 0.093
Adjusted R²: 0.0922
RMSE: 1.968
MSE: 3.872
MAE: 0.888

Nun fügen wir noch die Variable “bmi” zum Modell hinzu

Code

y_label3 = 'length_of_stay'

features3 = [#'recording_duration',
            #'age',
            'asa',
            'bmi',
            #'op_duration',
            #'estimated_blood_loss', 
            #'sex_F',
            ]

X3 = df[features3]
y3 = df[y_label3]

X_train3, X_test3, y_train3, y_test3 = train_test_split(X3, y3, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train3, y_train3)
y_pred3 = reg.predict(X_test3)

print("R²:",r2_score(y_test3, y_pred3).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test3, y_pred3))*((len(X_test3)-1)/(len(X_test3)-len(X_test3.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test3, y_pred3, squared=False).round(3))
print("MSE:",mean_squared_error(y_test3, y_pred3).round(3))
print("MAE:",mean_absolute_error(y_test3, y_pred3).round(3))

R²: 0.092
Adjusted R²: 0.0909
RMSE: 1.969
MSE: 3.876
MAE: 0.872

Da sich unser adjusted R² verbessert hat fügen wir nun noch die Variable “age” zum Modell hinzu

Code

y_label4 = 'length_of_stay'

features4 = [#'recording_duration',
            'age',
            'asa',
            'bmi',
            #'op_duration',
            #'estimated_blood_loss', 
            #'sex_F',
            ]

X4 = df[features4]
y4 = df[y_label4]

X_train4, X_test4, y_train4, y_test4 = train_test_split(X4, y4, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train4, y_train4)
y_pred4 = reg.predict(X_test4)

print("R²:",r2_score(y_test4, y_pred4).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test4, y_pred4))*((len(X_test4)-1)/(len(X_test4)-len(X_test4.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test4, y_pred4, squared=False).round(3))
print("MSE:",mean_squared_error(y_test4, y_pred4).round(3))
print("MAE:",mean_absolute_error(y_test4, y_pred4).round(3))

R²: 0.105
Adjusted R²: 0.1037
RMSE: 1.954
MSE: 3.819
MAE: 0.871

Unser Adjusted R² verbessert sich weiter, sodass wir die Variable “sex_F” zum Modell hinzufügen

Code

y_label5 = 'length_of_stay'

features5 = [#'recording_duration',
            'age',
            'asa',
            'bmi',
            #'op_duration',
            #'estimated_blood_loss', 
            'sex_F',
            ]
X5 = df[features5]
y5 = df[y_label5]

X_train5, X_test5, y_train5, y_test5 = train_test_split(X5, y5, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train5, y_train5)
y_pred5 = reg.predict(X_test5)

print("R²:",r2_score(y_test5, y_pred5).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test5, y_pred5))*((len(X_test5)-1)/(len(X_test5)-len(X_test5.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test5, y_pred5, squared=False).round(3))
print("MSE:",mean_squared_error(y_test5, y_pred5).round(3))
print("MAE:",mean_absolute_error(y_test5, y_pred5).round(3))

R²: 0.105
Adjusted R²: 0.1032
RMSE: 1.954
MSE: 3.819
MAE: 0.871

Das Adjusted R² hat sich nicht weiter verbessert, wir liegen jedoch weiter unter dem Wert für das Modell mit allen Prediktoren, daher fügen wir die Variable “estimated_blood_loss” zum Modell hinzu und nehmen die Variable “sex_F” wieder heraus.

Code

y_label6 = 'length_of_stay'

features6 = [#'recording_duration',
            'age',
            'asa',
            'bmi',
            #'op_duration',
            'estimated_blood_loss', 
            #'sex_F'
            ]

X6 = df[features6]
y6 = df[y_label6]

X_train6, X_test6, y_train6, y_test6 = train_test_split(X6, y6, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train6, y_train6)
y_pred6 = reg.predict(X_test6)

print("R²:",r2_score(y_test6, y_pred6).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test6, y_pred6))*((len(X_test6)-1)/(len(X_test6)-len(X_test6.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test6, y_pred6, squared=False).round(3))
print("MSE:",mean_squared_error(y_test6, y_pred6).round(3))
print("MAE:",mean_absolute_error(y_test6, y_pred6).round(3))

R²: 0.091
Adjusted R²: 0.0888
RMSE: 1.97
MSE: 3.881
MAE: 0.82

Das Modell verbessert sich nicht, wir prüfen, ob der Austausch der Variablen “estimated_blood_loss” und “op_duration” das Modell verbessert.

Code

y_label7 = 'length_of_stay'

features7 = [#'recording_duration',
            'age',
            'asa',
            'bmi',
            'op_duration',
            #'estimated_blood_loss',
            #'sex_F'
            ]

X7 = df[features7]
y7 = df[y_label7]

X_train7, X_test7, y_train7, y_test7 = train_test_split(X7, y7, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train7, y_train7)
y_pred7 = reg.predict(X_test7)

print("R²:",r2_score(y_test7, y_pred7).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test7, y_pred7))*((len(X_test7)-1)/(len(X_test7)-len(X_test7.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test7, y_pred7, squared=False).round(3))
print("MSE:",mean_squared_error(y_test7, y_pred7).round(3))
print("MAE:",mean_absolute_error(y_test7, y_pred7).round(3))

R²: 0.144
Adjusted R²: 0.1419
RMSE: 1.912
MSE: 3.655
MAE: 0.833

Das Modell hat sich weiter verbessert. Wir nehmen nun noch die Variable “recording_duration” hinzu.

Code

y_label8 = 'length_of_stay'

features8 = ['recording_duration',
            'age',
            'asa',
            'bmi',
            'op_duration',
            #'estimated_blood_loss',
            #'sex_F'
            ]

X8 = df[features8]
y8 = df[y_label8]

X_train8, X_test8, y_train8, y_test8 = train_test_split(X8, y8, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train8, y_train8)
y_pred8 = reg.predict(X_test8)

print("R²:",r2_score(y_test8, y_pred8).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test8, y_pred8))*((len(X_test8)-1)/(len(X_test8)-len(X_test8.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test8, y_pred8, squared=False).round(3))
print("MSE:",mean_squared_error(y_test8, y_pred8).round(3))
print("MAE:",mean_absolute_error(y_test8, y_pred8).round(3))

R²: 0.143
Adjusted R²: 0.1407
RMSE: 1.913
MSE: 3.658
MAE: 0.835

Das Modell verschlechtert sich. Wir prüfen nun, ob es unser Modell verbessert, wenn wir die “op_duration” im Vergleich zum eben berechneten Modell hinaus nehmen, Hintergrund ist, dass die “op_duration” und “recording_duration” eine starke Interkorrelation gezeigt hatten.

Code

y_label9 = 'length_of_stay'

features9 = ['recording_duration',
            'age',
            'asa',
            'bmi',
            #'op_duration',
            #'estimated_blood_loss',
            #'sex_F'
            ]

X9 = df[features9]
y9 = df[y_label9]

X_train9, X_test9, y_train9, y_test9 = train_test_split(X9, y9, 
                                                    test_size=0.3,
                                                    random_state=42)

reg.fit(X_train9, y_train9)
y_pred9 = reg.predict(X_test9)

print("R²:",r2_score(y_test9, y_pred9).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test9, y_pred9))*((len(X_test9)-1)/(len(X_test9)-len(X_test9.columns)-1))).round(4))
print("RMSE:",mean_squared_error(y_test9, y_pred9, squared=False).round(3))
print("MSE:",mean_squared_error(y_test9, y_pred9).round(3))
print("MAE:",mean_absolute_error(y_test9, y_pred9).round(3))

R²: 0.145
Adjusted R²: 0.143
RMSE: 1.91
MSE: 3.65
MAE: 0.826

Wir erhalten also mit den Variablen “recording_duration”, “age”, “asa” und “bmi” die besten Ergebnisse. Und speichern dieses Modell daher unter dem Namen LinearesRegressionsmodell.

Code

joblib.dump(reg, "C:\Studium\Project-applied-statistics\models\LinearesRegressionsmodell_2023-01-16.pkl")

['C:\\Studium\\Project-applied-statistics\\models\\LinearesRegressionsmodell_2023-01-16.pkl']

Lasso-Regressionsmodell

Wir werden nun ein Lasso-Regressionsmodell mit unseren Daten erstellen. Zuerst werden wir daher die Daten zu standardisieren, sodass unser Lasso-Regressionsmodell bestmöglich performt. Bedeutet, wir entfernen den Mittelwert und skalieren die Variablen auf eine Einheitsvarianz. Dies führen wir für Trainings- und Testdaten durch. Im Anschluss definieren wir unser Modell und integrieren gleichzeitig eine Kreuzvalidierung.

Code

scaler = StandardScaler().fit(X_train[features]) 

X_train[features] = scaler.transform(X_train[features])
X_test[features] = scaler.transform(X_test[features])
reg = LassoCV(cv=5, random_state=0)
reg.fit(X_train, y_train)

LassoCV(cv=5, random_state=0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Wir berechnen nun den Alpha-Wert für unser Lasso-Regressionsmodell

Code

reg.alpha_

0.0007694571302528419

Code

# Fit the model to the complete training data
regla = Lasso(alpha=reg.alpha_)
regla.fit(X_train, y_train)

Lasso(alpha=0.0007694571302528419)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Wir lassen uns nun die Regressionsgerade unsereres Lasso-Regressionsmodells berechnen.

Code

# intercept
intercept = pd.DataFrame({
    "Name": ["Intercept"],
    "Coefficient":[reg.intercept_]}
    )

# make a slope table
slope = pd.DataFrame({
    "Name": features,
    "Coefficient": reg.coef_}
)

# combine estimates of intercept and slopes
table = pd.concat([intercept, slope], ignore_index=True, sort=False)

round(table, 3)

	Name	Coefficient
0	Intercept	0.531
1	recording_duration	0.074
2	age	-0.232
3	asa	0.657
4	bmi	-0.125
5	op_duration	0.192
6	estimated_blood_loss	0.583
7	sex_F	0.064

Wir berechnen nun die Gütemaße unseres Modells.

Code

y_pred = reg.predict(X_test)

print("R²:",r2_score(y_test, y_pred).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test, y_pred))*((len(X_test)-1)/(len(X_test)-len(X_test.columns)-1))).round(3))
print("RMSE:",mean_squared_error(y_test, y_pred, squared=False).round(3))
print("MSE:",mean_squared_error(y_test, y_pred).round(3)) 
print("MAE:",mean_absolute_error(y_test, y_pred).round(3))

R²: 0.115
Adjusted R²: 0.112
RMSE: 1.943
MSE: 3.777
MAE: 0.802

Code

joblib.dump(reg, "C:\Studium\Project-applied-statistics\models\LassoRegressionsmodell_2023-01-16.pkl")

['C:\\Studium\\Project-applied-statistics\\models\\LassoRegressionsmodell_2023-01-16.pkl']

Gradient boosting

Wir versuchen mit Hilfe von Gradient boosting unser Lineares Regressionsmodell zu verbessern.

Code

y_label10 = 'length_of_stay'

features10 = ['recording_duration',
            'age',
            'asa',
            'bmi',
            #'op_duration',
            #'estimated_blood_loss',
            #'sex_F'
            ]

X10 = df[features10]
y10 = df[y_label10]

X_train10, X_test10, y_train10, y_test10 = train_test_split(X10, y10, 
                                                    test_size=0.3,
                                                    random_state=42)

params = {
    "n_estimators": 500,
    "max_depth": 4,
    "min_samples_split": 5,
    "learning_rate": 0.01,
    "loss": "squared_error",
}

reg = ensemble.GradientBoostingRegressor(**params)
reg.fit(X_train10, y_train10)


print("R²:",r2_score(y_test10, reg.predict(X_test10)).round(3))
print("Adjusted R²:",(1-(1-r2_score(y_test10, reg.predict(X_test10)))*((len(X_test10)-1)/(len(X_test10)-len(X_test10.columns)-1))).round(3))
print("RMSE:",mean_squared_error(y_test10, reg.predict(X_test10), squared=False).round(3))
print("MSE:",mean_squared_error(y_test10, reg.predict(X_test10)).round(3)) 
print("MAE:",mean_absolute_error(y_test10, reg.predict(X_test10)).round(3))

R²: 0.404
Adjusted R²: 0.403
RMSE: 1.595
MSE: 2.544
MAE: 0.549

Die Werte haben sich weiter verbessert. Wir werden diese nun aber im nächsten Abschnitt noch einmal gegenüberstellen.

Code

joblib.dump(reg, "C:\Studium\Project-applied-statistics\models\Gradient-Booosting-LinearesRegressionsmodell_2023-01-16.pkl")

['C:\\Studium\\Project-applied-statistics\\models\\Gradient-Booosting-LinearesRegressionsmodell_2023-01-16.pkl']

Results

Durch die Datenanalyse konnten wir erkennen, dass unsere Vorhersagevariable “length_of_stay” eine sehr große Streuung besitzt. Zudem konnten wir Zusammenhänge der Variablen betrachten und unser Lineares Regressionsmodell dementsprechend bestmöglich aufbauen. Störfaktoren in Form von Korrelation zwischen den Prädiktorvariablen wie “sex_M” und “sex_F” (r = -1,0) und “op_duration” und “recording_duration” (r = 0,989) wurden aufgezeigt und konnten bei der Modellerstellung berücksichtigt werden. Über die Methode der Vorwärtseliminierung konnten wir erkennen, dass ein Modell mit den Variablen “recording duration”, “age”, “asa” und “bmi” die besten Ergebnisse liefert. Nun können wir uns noch einmal die Ergebnisse, in Form der berechneten Gütemaße, des besten Linearen Regressionsmodells, des Lasso Regressionsmodells und des gradient boosting Modells gegenüber ansehen.

Code

results = pd.DataFrame({'': ['R²','Adjusted R²','RMSE','MSE','MAE'],
                        'Lineare Regression': [(r2_score(y_test9, y_pred9).round(3)),
                        ((1-(1-r2_score(y_test9, y_pred9))*((len(X_test9)-1)/(len(X_test9)-len(X_test9.columns)-1))).round(4)),
                        (mean_squared_error(y_test9, y_pred9, squared=False).round(3)),
                        (mean_squared_error(y_test9, y_pred9).round(3)),
                        (mean_absolute_error(y_test9, y_pred9).round(3))], 
                        'Lasso-Regression': [(r2_score(y_test, y_pred).round(3)),
                        ((1-(1-r2_score(y_test, y_pred))*((len(X_test)-1)/(len(X_test)-len(X_test.columns)-1))).round(4)),
                        (mean_squared_error(y_test, y_pred, squared=False).round(3)),
                        (mean_squared_error(y_test, y_pred).round(3)),
                        (mean_absolute_error(y_test, y_pred).round(3))], 
                        'Gradient boosting': [(r2_score(y_test10, reg.predict(X_test10)).round(3)),
                        ((1-(1-r2_score(y_test10, reg.predict(X_test10)))*((len(X_test10)-1)/(len(X_test10)-len(X_test10.columns)-1))).round(3)),
                        (mean_squared_error(y_test10, reg.predict(X_test10), squared=False).round(3)),
                        (mean_squared_error(y_test10, reg.predict(X_test10)).round(3)),
                        (mean_absolute_error(y_test10, reg.predict(X_test10)).round(3))]      
                        })

results

		Lineare Regression	Lasso-Regression	Gradient boosting
0	R²	0.145	0.1150	0.404
1	Adjusted R²	0.143	0.1118	0.403
2	RMSE	1.910	1.9430	1.595
3	MSE	3.650	3.7770	2.544
4	MAE	0.826	0.8020	0.549

Bei der Gegenüberstellung der Modelle fällt auf, dass das Lineare Regressionsmodell anhand aller Gütemaße mit Ausnahme des Mean-Absolut-Error besser abschneidet als das Lasso-Regressionsmodell. Das Modell mit Gradient boosting fällt bei allen Gütemaßen am besten aus und ist daher zu bevorzugen.

Lineares Regressionsmodell:

Unser Modell sagt mit den Variablen “recording_duration”, “age”, “asa” und “bmi” die Aufenthaltsdauer von Patienten am besten voraus. Dennoch sehen wir in unserem Modell nur ein kleines adj. R², was für eine geringe Güte des Modells spricht. Nur etwa 14,3% der Variabilität der Aufenthaltsdauer werden durch das Modell erklärt.

Lasso Regressionsmodell:

Wir sehen ein kleines Adjusted R² unseres Modells, was eine geringe Güte des Modells bedeutet. Die Güte ist sogar noch etwas unter der des linearem Regressionsmodells. Etwa 11,2% der Variabilität der Aufenthaltsdauer wird durch das Modell erklärt.

Gradient Boosting:

Wir sehen weiterhin ein kleines Adjusted R² unseres Modells, was für eine geringe Güte des Modells spricht. 40,3% der Variabilität der Aufenthaltsdauer werden durch das Modell erklärt.

Das Modell mit Gradient boosting schneidet bei allen Gütemaßen besser ab, als die beiden anderen Modelle und ist daher diesen gegenüber zu bevorzugen.

Anhand des MAE, Mean Absolute Error, erkennen wir den Durchschnitt der absoluten Differenz zwischen den von unserem Modell prognostizierten Werten und den tatsächlichen Werten. In unserem besten Modell bedeutet dies also, dass unser Modell die durchschnittlich die Aufenthaltsdauer um 0,55 Tage überschätzt. Im Falle, dass unser Modell im Krankenhaus verwendet würde, würde man tendenziell die Dauer 0,55 Tage länger planen, als diese tatsächlich ist. Der MSE, Mean Squared Error, zeigt uns den erwarteten, quadratische Abstand des prognostizierten Wertes zum tatsächlichen. Bei unserem besten Modell sind das 2,5 Tage. Der RMSE, Root Mean Squared Error ist eine Erweiterung von MSE und ist definiert als die Quadratwurzel des mittleren quadratischen Fehlers. Diese Kennzahl erreicht bei unserem besten Modell einen Wert von 1,6 Tagen. Bedeutet unser bestes Modell liegt durchschnittlich 1,6 Tage von der tatsächlichen Aufenthaltsdauer entfernt.

Discussion + Conclusion

Anhand der im “Results”-Teil beschriebenen geringen Güte unseres Modells, des RMSE und MAE würden wir einem Krankenhaus nicht dazu raten mit diesen Variablen und diesem Modell zu versuchen die Aufenthaltsdauer der Patienten zu prognostizieren. Bei einer durchschnittlichen Aufenthaltsdauer von 0,55 Tagen (alle Daten) bzw. 0,51 Tagen (Nachdem wir Aufenthaltstage >50 heraus genommen haben) einen RMSE von 1,6 Tagen zu haben scheint uns definitiv nicht gut genug.

Mögliche Verbesserungen wären: Über Stratified sampling Gruppen bilden, um so Stichproben auszuwählen. Gruppen könnten beispielsweise sein:

• „Abteilungen“ in welchen die Patienten behandelt werden
• Art der Operation
• handelt es sich um eine geplante oder eine Notoperation

Weiterhin wäre es auch möglich gewesen Variablen mit einzubeziehen, welche über die Forschungsfrage hinaus gehen. Hier wären ebenfalls die Variablen „Abteilungen“ in welchen die Patienten behandelt werden, Art der Operation und Notoperation denkbare Variablen gewesen.